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A large number of flows with distinctive patterns have been observed in experiments and simula- 
tions of Rayleigh-Benard convection in a water- filled cylinder whose radius is twice the height. We 
have adapted a time-dependent pseudospectral code, first, to carry out Newton's method and branch 
continuation and, second, to carry out the exponential power method and Arnoldi iteration to cal- 
culate leading eigenpairs and determine the stability of the steady states. The resulting bifurcation 
diagram represents a compromise between the tendency in the bulk towards parallel rolls, and the 
requirement imposed by the boundary conditions that primary bifurcations be towards states whose 
azimuthal dependence is trigonometric. The diagram contains 17 branches of stable and unstable 
steady states. These can be classified geometrically as roll states containing two, three, and four 
rolls; axisymmetric patterns with one or two tori; three-fold symmetric patterns called mercedes, 
mitubishi, marigold and cloverleaf; trigonometric patterns called dipole and pizza; and less symmet- 
ric patterns called CO and asymmetric three-rolls. The convective branches are connected to the 
conductive state and to each other by 16 primary and secondary pitchfork bifurcations and turning 
points. In order to better understand this complicated bifurcation diagram, we have partitioned it 
according to azimuthal symmetry. We have been able to determine the bifurcation-theoretic origin 
from the conductive state of all the branches observed at high Rayleigh number. 

PACS numbers: 47.20.Ky, 47.20.Bp, 47.10. Fg, 47.11. Kb 



I. INTRODUCTION 



In the late 1990s, Hof, Lucas and Mullin pj, |2| described five distinct steady patterns observed experimentally in 
a cylindrical Rayleigh-Benard convection cell at identical parameter values. More precisely, the patterns observed 
were torus, two-, three-, and four-roll states, and a mercedes pattern, at Prandtl number Pr = 6.7, Rayleigh number 
Ra = 14 200, and an aspect ratio F =radius/height=2 with insulating lateral boundaries. In our previous work [3,0], 
we reproduced numerically the five patterns observed by Hof and determined the approximate limits in Rayleigh 
number over which they could be observed. At lower Rayleigh numbers, we simulated several other patterns - dipole, 
pizza, and two-tori - as well as some time-periodic patterns. These results are summarized in figure [H Our viewpoint, 
pioneered in the 1980s by Benjamin and Mullin [5], is that these observations can be best understood and organized by 
constructing the bifurcation diagram corresponding to this figure. In particular, we wish to trace connections between 
the patterns observed at high and at low Rayleigh numbers, and to the basic conductive state wherever possible. 

The classical analysis of onset of Rayleigh-Benard convection describes an instability of the conductive state to 
a pattern of straight parallel rolls of infinite length. However, such a pattern is clearly not realizable in a small- 
aspect-ratio cylinder. Rolls must be curved to fit into the container, as shown in the two-, three- and four-roll states 
illustrated in figure [H In addition, a primary bifurcation, that is, a bifurcation from the conductive state, is associated 
with an eigenmode which is necessarily trigonometric in the azimuthal angle, such as the dipole or pizza states, or 
two-tori and torus states of figure [H The focus of this paper is the relationship between trigonometric modes and roll 
states and, more generally, the bifurcation-theoretic genesis of the profusion of states in this configuration. 

In our companion paper [4], we reviewed some of the literature on Rayleigh-Benard convection in small-aspect-ratio 
cylindrical geometries, focusing on pattern competition. The previous investigations most relevant to this manuscript, 
in addition to those of Hof et al. [l|, Q, are the full nonlinear simulations of Leong [6] and of Ma et al. [7]; we will 
compare our results to these articles where appropriate. 

In section |TT] we state the governing equations and the symmetries of the configuration. Section [111] describes the 
numerical methods we have used to compute steady states and their stability. Section IIVI begins by presenting the 
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full bifurcation diagram and primary bifurcations. We then give a detailed analysis of branches corresponding to each 
azimuthal wavenumber. Concluding remarks are presented in section IVl 




FIG. 1: (Color online) Schematic diagram of existence ranges and transitions between convective patterns observed in time- 
dependent simulation for insulating sidewalls. Stars denote solutions obtained from a slight perturbation of the conductive 
state, at the Rayleigh numbers indicated. The initial condition was identical for all five simulations. Arrows indicate patterns 
obtained by starting from stable steady states, and abruptly either lowering or raising the Rayleigh number. For example, at 
Ra = 2000, the perturbed conductive initial condition leads to a pizza state. Using the pizza state at Ra = 2000 as an initial 
condition leads to a four-roll state at Ra = 5000, to the three-roll state at Ra = 10 000, and to a two-roll state at Ra > 15 000. 
Right: representative patterns illustrated via temperature field in the horizontal midplane, with light portions representing hot 
rising fluid and dark portions representing cold descending fluid. 



II. BACKGROUND 



A. Governing equations and boundary conditions 

We recall from our companion paper [4] the dimensionless Navier-Stokes and Boussinesq equations governing the 
system: 

dtH^{\J-V)H = RaU.^V^H, (la) 
Pr-^(atU + (U- V)U) = -VP + V^U + i^e^, (lb) 
V-U = 0, (Ic) 

where H is the nondimensionalized deviation of the temperature from the linear vertical conductive profile. The 
parameter values are as follows: 

-p O 1 1 1 Q 

Pr = 6.7, r=-- — = 2, 0<i?a< 30000. (2) 

height ' - - \ J 
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The container is assumed to have rigid wahs, with thermahy conducting horizontal bounding plates and thermally 
insulating sidewalls 

U = for z = ±l/2 or r = T, (3a) 
H = for z = ±l/2, (3b) 
drH = for r = r (3c) 



B. Symmetries 

The bifurcations that this system can undergo are dictated by its symmetries. In group-theoretic terms, the 
conductive state has 0(2) symmetry in the azimuthal angle, meaning that it is invariant under all rotations and 
reflections in 0: 

{Ur,Ue,U,,H){r,e,z) = {Ur,Ug,U,,H){r,0 + eo,z) (4a) 
{Ur,Ue,U,,H){r,0,z) = {Ur, -Ue,U,, H){r,eo - 9, z) (4b) 

where Oq indicates an arbitrary angle of rotation or axis of reflection, and all compositions of these transformations. 
Under the Boussinesq approximation, the conductive state is also invariant under simultaneous reflection in z and 
change in sign of the temperature perturbation: 

{Ur, Ug, U,, H){r, e, z) = {Ur, Ue, -U,, -H){r, 9, -z) 

This symmetry can be combined with the ^-rotation symmetry (|^a|) to yield: 

{Ur, Ue, U„ H){r, 9, z) = {Ur, Ue, -U,, -H){r, 9 + 9o, -z) (4c) 

a form whose utility will appear shortly. The full symmetry group of the conductive state is thus 0(2) x Z2. 

A steady bifurcation from the axisymmetric conductive state, i.e. a primary bifurcation, is necessarily associated 
with an eigenvector which is trigonometric in the azimuthal direction; see, e.g. Crawford & Knobloch [14]. Each 
bifurcating branch is thus associated with an azimuthal wave number m. For m = 0, symmetry (|4c|) is broken and the 
bifurcation is a pitchfork, leading to two branches. If m is non-zero, the bifurcation is a circle pitchfork, producing 
families of states of arbitrary orientation. For the bifurcating states, 0(2) symmetry is replaced by Dm^ meaning that 
they are invariant under rotation by angles which are multiples of 27r/m and reflections in 2m axes of symmetry: 

{Ur,Ue,U,,H){r,9,z) = {Ur,Ue,U,,H){r,9 + 2n/m,z) (5a) 
{Ur,Ue,U,,H){r,9,z) = {Ur, -Ue,U,, H){r, -9, z) (5b) 
{Ur,Ue,U,,H){r,9,z) = {Ur,Ue,-U,,-H){r,9 + n/m,-z) (5c) 

where = is taken to be one of the axes of symmetry of the pattern, and ([5a|) is trivially verified if m = 1. These 
equations generate the symmetry group Dm x Z2. These states have a zero eigenvalue, corresponding to the marginal 
stability to rotation of the pattern. Equations (j5j) can be seen to be special cases of (The form of ([5c|) is the 
reason we choose (|4c|). instead of the Boussinesq reflection operator, as a generator of the symmetry group.) 

Primary branches can undergo secondary pitchfork bifurcations which break the Z2 symmetry ([5c|) . The result- 
ing branches, which we will call "asymmetric", nonetheless have Dm symmetry, generated by the discrete rotation 
symmetry ([5a|) and the reflection symmetry ([5b|) . 



III. NUMERICAL METHODS 



In [3|, we described our code for integrating the time-dependent Boussinesq equations in a cylindrical geometry. 
We have modified this time-dependent code using the techniques described in |8|, 19[ to carry out continuation by 
Newton's method and linear stability analysis by the exponential Arnoldi method. We describe these modifications 
in the subsections which follow. To do so, we will write the Boussinesq equations schematically as 

^=CU^^{U) (6) 

where C represents the viscous, diffusive and buoyancy operators and the advective terms. U = (H^Ur^Ue^Uz) 
represents the spatially discretized temperature deviation H and velocity field U = {Ur^Ue^Uz)- The imposition of 
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boundary conditions and incompressibility are assumed to be included in the representations of £, M and U. Here, 
we assume that timestepping is carried out via the first-order formula: 

U(t + At) = (/ - AtC)-\l + AtJ\f) U(t) = BiU{t)) (7) 

i.e. the terms in C are treated via the implicit backwards Euler scheme and those in M by the explicit forwards Euler 
scheme. 



A. Spatial discretization 

The code uses a pseudo-spectral spatial discretization, in which Uz are approximated as: 

Ne/2 2Nr--l 

f{r,e,z)=J2 E E W^j(^/r)rfc(2^)e™^ + c.c. (8) 

m=0 j>m. k=0 

j-\-m even 

while j+m odd is used for /7r, Uq. Differentiation is carried out on the spectral representation (|8]), while multiplications 
are performed after transforming to a grid, and then transforming the result back to the spectral representation. 
For the aspect ratio F = 2 investigated here, we use Nr = 40 gridpoints or Chebyshev polynomials in the radial 
direction, Nq = 120 gridpoints or trigonometric functions in the azimuthal direction and A^^ = 20 gridpoints or 
Chebyshev polynomials in the axial direction. Thus the domain is represented by approximately 10^ gridpoints and 
each solution by a vector of size 4 x 10^. (We have also checked our resolution for Ra > 20 000 by re-calculating 
a few of our branches - the mercedes, one-torus, two-roll and asymmetric three-roll branches - with a resolution of 
A/"^ X A^5> X A/'^ = 60 X 160 x 30.) The boundary conditions are imposed via the tau method, and incompressibility to 
machine accuracy is insured via an influence matrix technique. 



B. Steady state solving 

Steady states are found by calculating the roots of S — /, which are the same as those of A/" + £ for any value of 
At, as shown by the following calculation: 

(B-I) = {I - AtC)-\l ^ AW) - I 

= {I- AtC)-^ [(/ + AtN) -{I- AtC)] 

= {I-AtC)-^At{M^C). (9) 

The roots of B — I are found by Newton iteration: 

{Bu-I)u = {B-I)U (10a) 
U ^ U-u, (10b) 

where the Hnear operator Bu — I '^^ the Jacobian of B — I evaluated at U: 

[Bu - I)u = {I- AtC)-^At{Mu + C)u (11) 

while U = (i^, U) is the current estimate for the steady state and u = (/i, u) is an unknown correction to U. The 
action Muu is obtained from M{U) merely by carrying out the replacements 

VVH V -Vh^u-VH (12a) 
U-VU ^ U-Vu + u-VU (12b) 

in the nonlinear terms of ([Tb|) - ([Ta)) . Since the boundary conditions (j3j) are homogeneous, they remain unchanged. 

We iterate (p!Q]) until \ \{B — I)U\ \ is lower than some threshold, which we usually take to be CNewton = 10~^^, or until 
a maximum number of iterations, which we take to be 10, has been surpassed, meaning that the Newton procedure 
has failed. We use the norm 

^ ^max(||C/,||oo,||C/e||oo,||C/.||~))V (13) 



RaAt V Pr 
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The size of the matrix representing the hnear operator in (jlOap is (4 x 10^) x (4 x 10^) and so the system is far 
too large to be solved directly. Instead we use the BiConjugate Gradient Stabilized algorithm [10], which requires the 
right-hand-side and a procedure for calculating the action of — / on a vector u. The right-hand-side of (jlQap is 
shown by (|9|) to be the difference between U{t-\-At) = B{U{t)) and U{t), i.e. between two (widely spaced) consecutive 
timesteps, while the left-hand-side is the difference between Bu{U{t)) and U{t), i.e. between two /mean^^et/ timesteps. 
Conjugate gradient iteration proceeds until 

\\{Bu-I)u-{B-I)U\\ ^ 

— m-Duw — - '^'^^^ 

where the threshold CBiCGS is taken between 10~^ and 10~^^. The reason for finding the roots of S — / instead of 
those of A/" + £ is that, as shown by equation (|9|), 

{Bu -I)- C-^ {Mu + C) for At > 1. (15) 

This effective preconditioning by makes Bu — I far better conditioned than Nu + £, and greatly accelerates the 
convergence of BiCGSTAB. Note that At 1 is the limit opposite to that used in timestepping. We use At ranging 
between 0.2 and 10 (in contrast to the At on the order of 10~^ used in temporal integration). 

It is the solution of the linear system (|10ap which poses the greatest difficulty and which determined the limits of 
our study. In some regions, convergence of BiCGSTAB required as few as 5 actions of Bu — /, with more typical 
values ranging between 30 and 800. In other regions, 4000 iterations did not suffice (even sometimes far from any 
bifurcation, where singularity of Bu — / is to be expected), and continuation of the branch was eventually abandoned. 



C. Branch following 



In order to calculate a branch of steady states, we carry out Newton iteration (p!Q|) repeatedly for different values 
of Rayleigh number. Generally, in the absence of turning points, one can merely use the converged solution for one 
Ra to initialize the Newton iteration for a neighboring Ra. This initialization procedure constitutes zero-th order 
extrapolation. We reduce the increment or decrement ARa in Ra if the Newton iteration failed to converge in A/'°p* 
iterations and increase ARa if convergence took place sooner. Specifically, if we have computed solutions V(^^\ U^'^^ 
corresponding to Ra^^\ Ra^'^^ in N^^\ N^'^^ Newton iterations, we set 



7V°P* - 



ARa = Ra^'^^ + a{Ra^'^^ - Ra^^^) 



(16) 



1 



iV(2) + 1 



where we take 7V°p* between 2 and 5. 

Linear or quadratic extrapolation in Ra is easy to implement. Assume that converg ed solutions U^^\ U^^\ ZY^^) 
have been found for Rayleigh numbers Ra^^\ Ra^^^ and Ra^'^\ We can determine coefficients a^, 6^, q such that 



Ui = aiRa + hiRa + q 



(17) 



where i ranges over both the gridpoints and the components (i^, Ur^Uo^Uz). We then use (fT7|) to compute an initial 
condition for Newton iteration at the new value Ra'^^^ given in ([16]). (The condition number of the 3x3 system (fT7|) 
for a^, bi^ Ci is improved if we subtract from Ra the average of the three Ra values.) Over many portions of many 
branches we find we can easily take ARa > 200. As an example, we computed the marigold branch which will be 
described in section HV El from Ra = 2100 to Ra = 18 000, with intervals ARa varied dynamically between 10 and 
1200 according to prescription ([16]), requiring a computation time of 1200 CPU seconds on the NEC SX-8. 

Quadratic extrapolation, while not necessary for moving along a branch, proves essential near a turning point. Near 
a turning point {Ra^^ ,U^^ , we stop using extrapolation in Ra, as in ([T6]) . and instead use extrapolation in one of 
the components of U. That is, we fix the value of one component, Uj, and treat Ra as a dependent variable. To 
determine whether we are near a turning point, we use the fact that 



TP\ 



Ra-Ra^P\ 



(18) 



so that AUi must eventually exceed ARa as a turning point is approached. We monitor the relative changes by 
comparing the quantities 



m 










up 



with 



ARa 




Ra 





RgW - J?a(2) 
i?a(2) 



(19) 
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where 7 is a multiphcative weighting factor ranging between 5 (to favor extrapolation in Ra) and 0.001 (to favor 
extrapolation in Ui). When \/SlAi/Ui\ exceeds 7|Ai?a/i?a|, we replace ([16]), prescribing extrapolation in Ra^ by 
extrapolation in Ui\ 

Uf^ = up + AUi = up + a{uP - uP). (20) 

We use the three previous converged fields and Rayleigh numbers to determine coefficients a^, 6^, q for z ^ / and aRa^ 
bRa, CRa such that 

Ui = aiUj ^ hi Ui ^Ci, Ra = aRa Uj + hRa Ui + CRa (21) 

and then use ([2T]) to compute a new Ra and Ui^i ^ I corresponding to the Ui prescribed by (|2Q]) . This allows us to 
change direction in Ra] (|2T]) may lead to Ra^^^ — Ra^'^^ of opposite sign to that of Ra^'^^ — Ra^^\ unlike in equation 

m- 

The procedure above treats Ui as an independent variable (as prescribed in equation (|2Q|) ) and Ra as a dependent 
variable (as prescribed in equation (|2T]) ) in the predictor step (initialization). In this investigation, we have left the 
corrector step (Newton iteration) unchanged, that is, Ra remains unaltered by (|10ap . One strategy we have employed 
is to relax the tolerances near the turning point, for example to CNewton — 

10 and e^cGS = 10 ^. Like Xin 
we have succeeded in traversing a number of turning points in this way, despite the near-singularity of the matrix 
{Bu — I) near a bifurcation point. 



D. Linear stability analysis 

Once branches have been computed, we wish to determine their stability. In order to perform linear stability 
analysis of a steady state U = (i^, U), we carry out time integration of the Boussinesq equations linearized about U 
for an infinitesimal perturbation u = (/i,u): 

^=Cu^JVuu (22) 

We use the same timestepping formula ([7]) as for the nonlinear problem: 

u{t^At) = {I - AtC)-\l ^ AtAfu)u{t) =Buu{t) (23) 
by carrying out the substitutions in (fT2]) . Since 

Bu - exp(At(/: + Mu)) for At < 1, (24) 

the eigenvalues A of Bu and eigenvalues a of £ + Mu are related via 

A ^ exp(crAt) cr ^ log |A| for At <Cl (25) 

The stability of U is determined by the sign of the leading eigenvalue (Jmax (that with largest real part) of >C+A/L/, which 
corresponds to the dominant eigenvalue Amax (that with largest magnitude) of Bu- Equation ([23]) prescribes acting 
with the linear operator Bu on w(t); when repeated over many At's, u will converge to the eigenvector corresponding 
to Amax, which is itself approximated by the Rayleigh quotient 

_ {u{t),Buu{t)) 

Amax ^ lim . 26) 

To determine several leading eigenvalues, the power method is generalized to the Arnoldi-Krylov method [ll|. This 
consists of orthonormalizing a small number of fields 

{u{t = 0), u{t = T), u{t = 2T) . . . u{t = {K- 1)T)} (27) 

to create vectors vi^V2^vs^ . . .vk , and then a small Hessenberg matrix Hjk = (vj^BuVk), which is directly diagonalised. 
Its eigenvalues approximate eigenvalues A of Bu, while its eigenvectors (j) consist of coefficients of the vectors Vj, to 
be combined to form approximate eigenvectors ^ = (j)jVj of Bu- The accuracy of these approximate eigenpairs is 
measured by the residue | IS^^^ — A<l>| | in the case of complex eigenvalues. The integration of (|23l) is continued until 
the residues of the desired eigenvalues are below some acceptance criterion, usually near 10~^. 
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The timestep required is similar to that for timestepping. One obvious restriction comes from the exphcit scheme 
used in (|23|) for Mu ; a timestep which violates this stability requirement leads to approximate eigenvalues of Bu which 
bear no resemblance to those of exp(At(£ Mu))- For smaller At, the accuracy of the eigenvalues depends on At 
because of the approximation In particular, the time-splitting error means that Bu is not a function of C -^Mu- 
(In contrast, the errors in the pitchfork and turning point bifurcation thresholds obtained by Newton's method result 
only from the spatial discretization.) We have used At = 10~^ for Ra < 10 000 and At = 5 x 10~^ for Ra > simlO 000. 
We estimate our accuracy in locating bifurcation points and stability ranges to be ARa < 1. We used i^T = 10 vectors 
and a time interval of T = 100 At, i.e. T = 0.1 or T = 0.05 to create the Krylov vectors (|27|) . and an acceptance 
criterion for the residues of 10~^. 

A method which produces approximate eigenvalues which are independent of At is the inverse Arnoldi method [9[ , 
in which (|23|) is replaced by 

u"+^ = (/:+A/L/)-^u". (28) 

This is accomplished in practice by solving the equation 

{Bu - I)u''^^ = At(/ - At/:)-^u^ (29) 

The equivalence between (|28l) and (|29|) follows from a calculation similar to (|9]) . Equation (|29|) is very similar to (jlOap 
and is also solved by BiCGSTAB. Only a few iterations (between 1 and 10) of (|29]) lead to an extremely accurate 
eigenvalue. However, the inverse Arnoldi is more difficult to implement than the exponential Arnoldi method. For 
this reason, we have chosen not to do so for this study. 



IV. RESULTS 



A. Bifurcation Diagram 

Using the methods described in section IIIIl we have succeeded in continuing the branches we found previously 
via time integration [sl, B]. By going around turning points and bifurcation points, we have computed a total of 17 
branches of convective steady states. These are related to the conductive state and to each other by 5 primary and 3 
secondary pitchfork bifurcations, and 8 saddle- node bifurcations. The bifurcation diagram is shown in figure [2j Tables 
summarizing all of the branches and bifurcations we have found are given in section |Vl 

The axes of figure [2] have been chosen with care. In order to show the full extent of our calculations in Rayleigh 
number and, at the same time, distinguish between various low- Rayleigh- number primary bifurcations, figure [2] uses a 
logarithmic scale in Ra. More specifically, using \og{Ra — 1000) distinguishes primary bifurcations better than log Ra. 
The vertical axis was chosen to best distinguish between the various branches. The quantity H is the maximum 
absolute value of the temperature deviation over the ring at (r = 0.3, O^z = 0) 

H = max|i^(r = 0.3,6>,z = 0)|. (30) 

9 

H itself and the commonly used Nusselt number deviation 

Nu-1 = J rdr dO d^H(r, 0,z = 0) - 1 (31) 

have a strong linear dependence on Ra; plotting them directly as a function of Ra does little to separate the branches. 
We have therefore chosen instead to represent each state by its value of H/Ra. (Exceptionally, for the first two-tori 
branch, we have plotted —H/Ra for low Ra to avoid a reversal in slope due to the absolute value in (|3Q|).) Each 
branch in figure [2] is a representative of a number of branches - the group orbit - that can be obtained by refiection 
and rotations. 

In order to understand the complicated bifurcation diagram in figure [2l we will select various aspects for detailed 
study below. 



B. Primary Bifurcations 



We give in table [J the first critical wavenumber and Rayleigh number pairs. The thresholds given to ARa = 0.1 
are extrapolations from the branches calculated using Newton's method. The thresholds given as integer values were 
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FIG. 2: (Color online) Bifurcation diagram containing 17 branches of steady states, in addition to the conductive branch 
(indicated by short-dashed horizontal line). Shown are pizza (solid green), four-roll (long-dashed turquoise), two-tori (solid red; 
2), torus (long-dashed magenta; 2), marigold (solid blue), mitsubishi (short-dashed purple), cloverleaf (long-dashed purple) 
and mercedes (solid blue), three-roll (solid black), tiger (long-dashed brick), asymmetric three-roll (solid brick; 2), two-roll 
(solid blue; 2), CO (long-dashed red) branches. The notation "torus (long-dashed magenta; 2)", e.g., signifies that there are 2 
torus branches, related by saddle-node bifurcations and both shown as long-dashed magenta curves. Note that the bifurcation 
diagram gives no information concerning stability; i.e. whether a solution curve is depicted as solid or dashed does not indicate 
its stability. Turning points or pitchfork bifurcations are shown as dots. 



calculated from the linear stability analysis of the conductive branch. In the remainder of the manuscript, we round 
Rayleigh numbers to integer values (except in a few very specific cases). Our thresholds agree quite well with those 
of previous researchers. The discrepancies are typically on the order of 0.3% with the calculations of Ma et al. 0] 
and on the order of 0.02% with those of Martin- Witkowski [13], which we believe to be the two most recent threshold 
calculations in this geometry. With increasing Rayleigh number, many other bifurcations occur from the conductive 
state, both to higher wavenumbers and to different eigenmodes with the same wavenumbers. The branches created 
at such bifurcations are necessarily unstable. 

It is the first four bifurcations of table [H along with the last column, which will prove relevant to the steady states 
observed, i.e. the stable ones. In figure [3] we show these first bifurcations, along with corresponding nonlinear states 
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TABLE I: (Color online) First bifurcations from conductive state 
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FIG. 3: (Color online) Primary branches bifurcating from conductive state. For this aspect ratio, F = 2, and with insulating 
lateral walls, the first four critical wavenumber and Rayleigh number pairs are (m = 1, Ra = 1828; black), (m = 2, Ra = 1849; 
green), (m = 0, Ra = 1861; red), and (m — 3, Ra — 1985; blue). Below are representative states from each of the bifurcating 
branches. 

at slightly supercritical values of Ra. We recognize the dipole, pizza and two-tori states. The other states in figures 
[Hand [2]- the two-, three-, and four-roll states, or the torus, mercedes, and CO states - are not present in figure [3l 
Their origin is addressed in the sections which follow. 



C. Pizza and Four- roll branches (m = 2) 

We now focus on various sets of solution branches. We begin with the branches arising from the instability to an 
m = 2 quadrupolar eigenvector at Ra = 1849, because these are free from the complications which we will encounter 
for the other azimuthal wavenumber s. We use three figures to describe the structure of these branches. Figure |4] 
uses the same coordinates as figure [2l merely extracting the relevant branches. Figure [5] is a qualitative bifurcation 
diagram, accompanied by illustrations of representative states along the branches. Finally, figure [6] shows leading 
eigenvalues, from which the stability of the underlying branches can be deduced. 

The bifurcation sequence is best understood by studying figure [5l The schematic quantity along the vertical axis 
and the monotonic but non-uniform Rayleigh- number progression along the horizontal axis, are chosen to separate 
the different branches and to illustrate the bifurcations. Representative states along the branches are illustrated 
via temperature distributions in the midplane {z = 0), with light portions representing hot rising fluid and dark 
portions representing cold deseeding fluid. To avoid further cluttering the figure, the Rayleigh numbers given for the 
representative states have been rounded to the nearest 10, 100 or even 1000, with precise bifurcation points given 
along the axis. The azimuthal orientation of the representative states is arbitrary, and to each branch corresponds 
another branch obtained by the Boussinesq reflection symmetry which, for these illustrations, would mean reversing 
light and dark. 
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FIG. 4: (Color online) Partial bifurcation diagram including m — 2 primary branch and connecting branches. Bifurcations are 
shown as dots. The primary Pizza branch (solid, green) bifurcates from the conductive state at Ra — 1849 and terminates 
in a saddle-node bifurcation at Ra ~ 19 450. The Four-roll branch (dashed, turquoise) bifurcates from the pizza branch at 
Ra = 2353 and terminates at a saddle-node bifurcation at Ra ^ 23 130. 
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FIG. 5: (Color online) Schematic partial bifurcation diagram relating branches originating from the m — 2 bifurcation. At 
Ra — 1849 the Pizza branch originates via a circle pitchfork bifurcation from the conductive state corresponding to an m = 2 
eigenvector. It terminates at a turning point at Ra < 19 450 and is stable for 1879 < Ra < 2353. At Ra = 2353, a secondary 
pitchfork bifurcation leads to a Four-roll branch, which is stable for Ra < 22 660 and ends at a turning point at Ra ~ 23 130. 
For visual clarity, the Rayleigh numbers given for the representative states have been rounded to the nearest 10, 100 or 1000. 



A circle pitchfork bifurcation from the Conductive branch to an m = 2 eigenmode takes place at Ra = 1849. 
Figure [5] shows that, near onset, the states along the branch created by this bifurcation contain two hot upwelling 
spots and two cold downwelling spots. Their resemblance to a small pizza leads us to call this the Pizza branch. 
As Ra increases, the central convective regions shrink. By the time the pizza branch terminates at a saddle-node 
bifurcation at Ra = 19 450, most of the convection takes place at four regions along the edge of the container. 

A pitchfork bifurcation at Ra = 2353 from the pizza branch breaks the symmetry between hot upwelling and 
cold downwelling fluid: the two downwelling spots merge as the two upwelling spots elongate (or vice versa for the 
complementary branch, not shown). This secondary bifurcation is also computed by Ma [7], who gave its threshold 
as Ra = 2350. The pitchfork bifurcation leads to a Four-roll branch which terminates at a saddle- node bifurcation 
at Ra ~ 23 130. Along the four-roll branch, the convective regions diminish as Ra increases, as was the case for the 
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FIG. 6: Four leading eigenvalues of (a) the pizza branch at low Ra and of (b) the four-roll branch at high Ra. Bifurcations 
(zero crossings) indicated by dots. The zero eigenvalue (dotted) which exists throughout corresponds to the marginal stability 
to rotation of the pattern, a) The bifurcating eigenvalue (short-dashed) decreases steeply from at onset, Ra = 1849. The 
pizza branch is initially unstable since it inherits the unstable eigenvalue (solid) of the conducting branch, due to the preceding 
m = 1 bifurcation. This leading eigenvalue decreases with Ra, crossing zero at Ra = 1879. Another eigenvalue (long-dashed) 
becomes positive at Ra = 2353, accompanying the bifurcation to the four-roll branch. The stability interval of the pizza branch 
is 1879 < Ra < 2353. b) The four-roll branch loses stability near Ra = 22 660. 



pizza branch; the rolls become wide, with narrow upwelling and downwelling boundaries. 

Figure [6]a shows the four leading eigenvalues of the pizza branch near onset, computed by the methods described 
in section IIIIDI They are grouped into distinct sets by examining the spatial structure, especially the azimuthal 
wavenumber spectrum, of the corresponding eigenvectors, and then plotted as curves. Very near onset, each eigen- 
value can be associated with an azimuthal wavenumber, since it is connected to an eigenvalue of the conductive 
branch. The zero eigenvalue (m = 2, sometimes called the phase mode) which exists throughout corresponds to the 
marginal stability to rotation of the pattern. The eigenvalue which is zero at onset and then rapidly decreases is that 
corresponding to the circle pitchfork which creates this branch (also m = 2, sometimes called the amplitude mode). 
The positive eigenvalue (m = 1) at onset results from the fact that the bifurcation to the dipole branch at Ra = 1828 
precedes the creation of the pizza branch. The pizza branch inherits this instability when it is created at Ra = 1849 
and becomes stable at Ra = 1879, when the leading eigenvalue becomes negative, as shown in figure [6]a. This confirms 
our time-dependent simulations [3, 4], summarized in figured] which shows the pizza branch as stable at Ra = 2000. 
However, it enjoys only a short i?a-interval of stability, as another eigenvalue (connected to the m = eigenvalue of 
the conductive branch) becomes positive at Ra = 2353, when the secondary pitchfork bifurcation creates the four-roll 
branch. 

Figure [6]?) shows that the four-roll branch remains stable until Ra < 22 660, a far wider Rayleigh-number interval 
than the pizza branch. Indeed, roll states are preferred by convective systems and are those generally observed 
in experiments and time-dependent simulations. In particular, a four-roll state was computed in time-dependent 
simulations [3, 4] for Ra between 5000 and 20 000 (see figure [T|) and is one of the five states observed experimentally 
by Hof et al [1] at Ra = 14 200. 

As explained in section |IIB| the four-roll states have symmetry group D2: they are invariant under rotation by tt 
and reflection in either of the symmetry axes, as stated in (|5a|) - (|5bp . The pizza states are also invariant under the 
additional symmetry given in (|5c|) (rotation by 7r/2, z — Uz —Uz^ H —H) so have the larger symmetry 
group X ^2, as is typical for the primary branches bifurcating from the conductive state. 



D. Torus and two-tori branches (m = 0) 



We now survey the axisymmetric branches. Figure [71 extracts the axisymmetric branches from the complete bifur- 
cation diagram of figure [2l There are two pairs of branches, i.e. a total of four branches of axisymmetric states. A 
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schematic bifurcation diagram showing representative states is given in figure [8] and leading eigenvalues are shown in 
figure [9l 

The Two- tori branches result from two pitchfork bifurcations from the conductive state at Ra = 1862 and Ra = 
2328. These two branches meet and annihilate at a turning point at Ra = 12 711. Most of the states along these 
branches contain two concentric toroidal convection rolls. The branch created at Ra = 1862, which we call the upper 
or stable two-tori branch, is the more stable of the two. In fact, it is unstable when it is first created, as shown in figure 
[H since the m = 1 and m = 2 bifurcations precede the m = bifurcation. For an axisymmetric convective state, the 
eigenvectors are each associated with a single azimuthal wavenumber m. The bifurcating eigenvalue, with m = 0, is 
at onset and rapidly decreases. One of the two leading eigenvalues becomes negative at Ra = 2116 (m = 2) and the 
second at Ra = 2300 (m = 1), stabilizing this two-tori branch. These stabilizing bifurcations were also computed by 
Ma et al 0, with thresholds 2113 and 2245, respectively. The upper two-tori branch remains stable until Ra = 5438, 
when the m = 1 eigenvalue becomes positive again. 

The One-torus branches emerge from a saddle-node bifurcation at Ra = 3076. The states on these branches all 
contain a single toroidal convection roll. We have not found any connection between these branches and any others, 
including the conductive branch. Both branches are initially unstable. The lower one-torus branch never stabilizes 
and we have been unable to calculate it past Ra = 17 857. The upper or stable one-torus branch is created with five 
positive eigenvalues. As Ra increases, these successively become negative, as shown in figure [9l The branch is stable 
for Ra > 4918 and exists until at least Ra = 29 940. It is clear that the axisymmetric state observed at Ra = 14200 
in experiment [1] and in time-dependent simulation [3, 4, 6, 7] must be on the stable one-torus branch, and not on 
the two-tori branch (which is unstable for Ra > 5438 and does not exist for Ra > 12 711) which bifurcates from the 
conductive state. 

Along the upper two-tori branch the inner roll dominates, while along the lower two-tori branch, the outer roll 
dominates. These states do not necessarily all contain two rolls. In particular, some states along the lower two-tori 
branch for Ra < 5000 seem to contain only one roll. These states bear a qualitative and quantitative resemblance to 
those on the upper one-torus branch. On figure [71 the upper one-torus and the lower two-tori branches are nearly 
tangent to one another over the interval 3000 < Ra < 3500, while at the turning point at Ra = 3076, the one-torus 
states bear a strong resemblance to the lower two-tori branch. 

The axisymmetric convective branches break the Boussinesq symmetry (|4c|) . while retaining the 0(2) azimuthal 
symmetry (|^ - (|4b|) . 

E. Mercedes, Cloverleaf, Mitsubishi and Marigold states (m = 3) 

The set of branches with three- fold symmetry are perhaps the most interesting, and certainly the most aesthetic. We 
have been able to trace the tortuous connection between the states obtained by time-integration (and hence necessarily 
stable) and the m = 3 primary branch (which, occurring after three other primary bifurcations, is necessarily unstable). 
Figure [To] extracts from figure [2] the branches with three- fold symmetry. Four branches are present (in addition to the 
conductive branch), connected by saddle- node and pitchfork bifurcations, shown as dots. Figure [TT] plots the leading 
eigenvalues of these four branches. The qualitative bifurcation diagram in figure [12] provides a clearer picture of the 
bifurcations, without the crossings present in figure [10] We recall that an identical set of branches, with hot and cold 
reversed (along with upwelling and downwelling) also exists, and that the azimuthal orientation is arbitrary. 

We begin by describing the states in figure [121 beginning from the stable state at Ra ^ 30 000, which Hof called 
Mercedes because of its resemblance to the logo of this automobile [2]. This state was observed in Hof's experiment [l| 
and in time-dependent simulations [3, 4], where it was computed for Ra between 5000 and 29 000, as shown in figure 
[H A mercedes state was also computed at Ra = 31 250 by Leong Q and at Ra = 14 200 by Ma et al Q. For high 
Ra, the hot upwelling and cold downwelling regions in the midplane are narrow, confined to three hot spots along 
the lateral boundary and a central cold Y-shaped region. With decreasing Ra, the upwelling and downwelling regions 
come to occupy an increasing portion of the midplane. By the turning point at Ra = 4634, the three hot spots have 
widened, becoming almost circular, and the center and extremities of the cold Y have widened into four triangles. 
Emerging from this turning point is what we have called the Cloverleaf branch. Following this branch towards 
increasing values of Ra, the three hot spots move inwards from the boundary and the cold Y-shaped region breaks, 
leaving four separate triangles. The hot spots and the central cold triangle become smaller, while the three remaining 
cold triangles narrow and cling to the lateral boundary. By the turning point at Ra = 18 762, the hot spots have 
merged into one central triangular region, and the cold regions form a ring occupying almost the entire circumference. 
As we follow the new branch with decreasing Ra, the points of the triangle expand and separate, forming oval petals 
or blades, while the exterior ring forms three exterior triangles. Midway along this branch, the states resemble the 
logo of the Mitsubishi automobile, and this is the name we have given to the branch. The hot uprising and cold 
downwelling regions become more similar as Ra decreases. At Ra = 4103, the Mitsubishi branch is seen to emanate 
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FIG. 7: (Color online) Partial bifurcation diagram including only axisymmetric states. The Two-tori branches (dashed, 
magenta) emerge from pitchfork bifurcations from the conductive branch at Ra — 1861.5 and Ra — 2328. They join and 
terminate at a turning point at Ra — 12 711. The One-torus branches (solid, red) emerge from a turning point at Ra — 3076 
and seem to be unconnected to the conductive state. 
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FIG. 8: (Color online) Schematic partial bifurcation diagram showing axisymmetric branches. The Two-tori branches are 
connected to the conductive branch via pitchfork bifurcations at Ra — 1862 and 2328, and to each other via a turning point 
at Ra — 12 711. The upper two-tori branch is stable for 2300 < Ra < 5438. The One-torus branches are connected to each 
other via a turning point at Ra = 3076. The upper one-torus branch is stable for Ra > 4918. Most states along the two- tori 
branches contain two concentric rolls, but states on the lower two-tori branch resemble those on the upper one-torus branch 
for Re < 3500. 
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FIG. 9: Five leading eigenvalues of upper axisymmetric branches, each labelled with its azimuthal wavenumber m. a) The 
upper two-tori branch has two positive eigenvalues at onset at Ra = 1862, which cross zero at Ra = 2116 and Ra = 2300 and 
again at Ra = 5438. It is stable for 2300 < Ra < 5438. b) The upper one-torus branch has five positive eigenvalues at onset 
at Ra = 3076. These cross zero at Ra =3330, 3438, 4408, 4582 and 4918, above which the branch is stable. 

in a pitchfork bifurcation from the Marigold branch, whose states have six equal petal-shaped regions. The marigold 
branch itself is generated at a circle pitchfork bifurcation from the conductive branch at Ra = 1985. 

The Mitsubishi, cloverleaf and Mercedes states have Ds symmetry, while the marigold states have the larger 
symmetry group Ds x Z2. 

The cloverleaf and Mitsubishi branches were obtained from the Mercedes branch by going around the turning points 
via the quadratic extrapolation described in section IIIII Additional effort is required to switch from the Mitsubishi 
to the marigold branch, since straigtforward continuation treats the pitchfork as it would a turning point. Because 
we can calculate eigenvectors, steady states and transient behavior, there are a number of ways in which a starting 
point on the marigold branch could be obtained: 

(i) Take a Mitsubishi state, for which the left and right-hand-sides of (|5c|) are not equal, and average the two 
expressions. 

(ii) Add a small amount of the m = 3 eigenvector to the conductive branch. 

(iii) Carry out time-integration for Ra > 1985, retaining only azimuthal modes which are multiples of three. 

We used method (iii), halting the integration after the marigold state was reached but before its instability was 
manifested. 

Although figure [11] shows leading eigenvalues corresponding to the states of figure [TOl it is of a different nature 
from the previous eigenvalue plots. Figures [6] and [9] showed one or more leading eigenvalues for states along a single 
branch. In contrast, figure [11] shows a single eigenvalue per state, but for states along the four different branches 
described above. Thus, between one and four eigenvalues are shown for a single Rayleigh number. All of the branches 
have at least one positive eigenvalue (and are thus unstable) except the Mercedes branch for Ra > 5503. When the 
marigold branch bifurcates from the conductive branch at Ra = 1985, it inherits three positive eigenvalues, the largest 
of which is the highest curve in figure [11] The Mitsubishi branch shares the spectrum of the marigold branch at the 
pitchfork bifurcation at Ra = 4103. As Ra increases, both branches become more unstable, but, for Ra > 12 500, 
the eigenvalues of the Mitsubishi branch eventually begin to decrease. The Mitsubishi branch is still unstable when 
it meets the cloverleaf branch at the turning point at Ra = 18 762. Following the cloverleaf branch with decreasing 
Ra, the leading eigenvalue decreases. When the cloverleaf and Mercedes branches meet at the turning point at 
Ra = 4634, the leading eigenvalue is still barely positive. Following the Mercedes branch with increasing Ra, the 
leading eigenvalue continues to decrease, becoming negative and stabilizing the branch for Ra > 5503. 
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FIG. 10: (Color online) Partial bifurcation diagram including only branches with three-fold symmetry. The Marigold branch 
(solid, blue) arises at a pitchfork bifurcation from the Conductive branch (short-dashed, black) at Ra — 1985. A secondary 
pitchfork bifurcation from the marigold branch at Ra = 4103 gives rise to the Mitsubishi branch (short-dashed, lighter purple). 
At a turning point at Ra — 18 762, it meets the Cloverleaf branch (long-dashed, darker purple). The Mercedes branch (solid, 
blue) originates at another turning point at Ra = 4634. 




FIG. 11: (Color online) Leading eigenvalue for each of the three- fold-symmetric branches. From highest to lowest: marigold 
(solid, blue Mitsubishi (short-dashed, light purple), cloverleaf (long-dashed, dark purple), Mercedes (solid, blue). Dots indicate 
bifurcations from conductive to marigold branch {Ra = 1985, a ^ 1.46), to Mitsubishi branch {Ra = 4103, a ^ 10), to 
cloverleaf branch, {Ra = 18 762, a ~ 27), to Mercedes branch {Ra — 4634, a « 1-3), and final stabilization of Mercedes branch 
{Ra = 5503, a = 0). 
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FIG. 12: (Color online) Schematic partial bifurcation diagram relating branches originating from the m — 3 bifurcation. The 
four branches of steady states with three-fold symmetry are the Mercedes, Cloverleaf, Mitsubishi, and Marigold branches. 
The Marigold branch is created by an m = 3 circle pitchfork bifurcation from the Conductive branch at Ra = 1985. It undergoes 
a pitchfork bifurcation at Ra = 4103, leading to the Mitsubishi branch. A turning point at Ra = 18 762 leads to the cloverleaf 
branch, and another turning point at Ra — 4634 to the Mercedes branch. Only the Mercedes branch is stable, for Ra > 5503, 
as indicated by the thick line. 



F. Dipole, Tiger and Three-roll branches (m =1) 



Although the m = 1 bifurcation has the lowest Rayleigh number threshold, Ra = 1828, we have postponed its 
discussion because of its odd behavior. Two branches bifurcate simultaneously, as shown in figures [13] and [Ml The 
states close to the bifurcation are of dipole form, as expected. A dipole state is observed at Ra = 2500 by Leong Q 
and by Ma et al [7]. Along one branch, additional spots appear on either side of the dipole, which grow as Ra is 
increased; we have given the name of Tiger to this branch. We have been unable to compute the tiger branch above 
Ra = 7936. The other branch is more conventional. The two parts of the dipole elongate along the dipole axis and 
patches of opposite sign appear and elongate near the boundary, leading eventually to a Three-roll structure. The 
three-roll branch exists at least until Ra = 30 000, where we have stopped our computations. 

Time-dependent simulations P, 0] between Ra = 20 000 and 25 000 show a transition to an Asymmetric Three- 
roll state, for which the rolls are slightly shifted. We have determined that an asymmetric branch emerges from the 
three-roll branch at Ra = 22 125 via a subcritical pitchfork bifurcation (the only subcritical bifurcation we have found 
in this system). The branch reverses direction and stabilizes at a saddle- node bifurcation at Ra = 21078. Three-roll 
and asymmetric three-roll states are observed at Ra — 12 500 and Ra — 25 000, respectively, by Leong [6]. 

The tiger and the three-roll branches share the same spatial symmetry, that is D\ x Z2, where D\ (equivalent to 
Z2) is generated by reflection in the axis perpendicular to the roll or dipole axis, as in (|5b]), and Z2 is generated by 
simultaneous rotation by tt and reflections z — Uz —Uz^ H — i7, as in (|5c|) . The asymmetric three-roll 
branch has symmetry Di. 

Figure [T5]a, containing leading eigenvalues of each branch near threshold, shows that the tiger branch becomes quite 
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FIG. 13: (Color online) Partial bifurcation diagram of branches created at the m — 1 primary bifurcation. Two branches of 
dipole states are created at Ra = 1828, and seem to be tangent to one another near the bifurcation. As Ra increases, the spatial 
forms along these branches evolve in divergent ways. One branch contains the three-roll states (solid, black) and the other 
(dashed, brick) contains the tiger states. An asymmetric three-roll branch (solid, brick) is created at a subcritical pitchfork 
bifurcation at Ra = 22 155 and reverses direction at a saddle- node bifurcation at Ra — 21 078. 
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FIG. 14: (Color online) Schematic bifurcation diagram of branches emanating from m — 1 instability. Two branches bifurcate 
simultaneously at Ra — 1828. Near the bifurcation, states along both branches have the form of a dipole. States along one 
branch evolve into a Three-roll state, while along the other branch evolve to a form called Tiger. The three-roll state is 
stable for 3762 < Ra < 20 393, and the tiger branch is never stable. The Asymmetric Three- Roll branch is formed at 
Ra — 22 125, via a subcritical pitchfork bifurcation from the three-roll branch, reversing direction at a saddle-node bifurcation 
at Ra = 21 078. 
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FIG. 15: (Color online) a) Leading eigenvalues of the tiger (long-dashed, brick) and three-roll (solid, black) branches. The 
eigenvalues of the tiger branch grow rapidly as Ra increases, while the three-roll branch is weakly unstable until Ra = 3762. 
The rapidly decreasing eigenvalue is that associated with the formation of these branches, b) The tiger and three-roll branches 
have the same curvature near onset. 



unstable immediately after it forms. At the same time, the three-roll branch also becomes weakly unstable, though 
it eventually stabilizes at Ra = 3762. Instability near onset is another unexpected feature of these branches. The 
three-roll branch becomes unstable again at Ra = 20 393, as shown in figure [T6l 

We have sought to better understand the primary m = 1 bifurcation, at which both the tiger and three-roll branches 
bifurcate simultaneously. First, we have verified that the three-roll and the tiger branches are distinct by following 
them around the pitchfork bifurcation, producing symmetrically-related branches. Thus the possibility that our 
continuation procedure has jumped from one branch to another is ruled out. Figure [T5b shows that both the tiger 
and the three-roll branches emerge via a pitchfork bifurcation, i.e. that H oc \/Ra^^^M^^ or, equivalently, that 

^^_\ _ — Rac . X 

Rac) ~^ Rac ^ 

for Rac = 1828.37 < Ra < 1900. Figure [T5b also shows that the constant of proportionality in (|32]) is the same 
for both branches {a ~ 0.59). Thus the two branches initially share not only a vertical tangent, but even the same 
curvature. 

Simultaneously bifurcating and non-equivalent branches are encountered in a number of situations, notably pitchfork 
bifurcation in the presence of D4 symmetry, such as in a square box [3, [HI. In the D4 case, one set of branches 
contains solutions whose axes of symmetry are the vertical or the horizontal midline of the square, while, for the other 
set of solutions, the symmetry axes are the diagonals. Although the two types of branches bifurcate simultaneously, 
they are not related to one another by a symmetry operation of D4, and so are not dynamically equivalent. For 0(2), 
in contrast to D4, the concepts of vertical, horizontal and diagonal have no meaning: solutions of any orientation can 
be obtained by rotation, and so must all be dynamically equivalent. An explanation of the m = 1 behavior is the 
subject of a separate investigation. 

Most of the results given above are consistent with the time-dependent simulations of our companion paper Q , 
summarized in figure [TJ The inconsistencies can largely explained as a consequence of finite integration time and very 
weak instability. This accounts for the observation of a symmetric three-roll state for 25 000 < Ra < 30 000, a range 
over which it is unstable, rather than its stable asymmetric counterpart. The same is true of the long-lived dipole 
state observed near onset at Ra = 2000. In all of these cases, the largest eigenvalue is less than 0.5, which would 
not have allowed initally small perturbations to grow to appreciable levels over the duration of our time- dependent 
simulations; this is true for all of the steady states shown in figure [T] summarizing the time-dependent simulations. 

Time-dependent simulation from an initial dipole state at various Rayleigh numbers also yielded various interesting 
transients [3|, and led to the discovery of an oscillatory pattern, the rotating S, and of two new steady patterns: 
the dipole smile, whose branch we have not continued, and the CO, described in the next section. 
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FIG. 16: (Color online) Enlargement of bifurcation diagram of figure [T3l The symmetric three-roll branch loses stability at 
Ra — 20 393. A subcritical bifurcation at Ra — 22 125 from the symmetric three-roll branch leads to the creation of the 
asymmetric three-roll branch, which is stabilized by a saddle-node bifurcation at Ra = 21 078. 



G. Two-roll and CO branches 

Finally, we mention another set of branches which, like the one-torus branches, appear to be unconnected to the 
conductive state. These are the Two-roll and the CO branch, shown in figures [171 and [TSl 

Figure [1] shows how these branches were originally found by time-integration, staring from quasi-steady states at 
Ra = 2000. A CO state was found by starting from a dipole and setting Ra — 10 000 and a two-roll state by starting 
from a pizza and setting Ra = 16 000. The two-roll branch originates at a saddle-node bifurcation at Ra = 8677, 
where it is connected to an unstable branch containing states which also have two rolls. Figure [18] shows that, for 
high i?a, states on the stable branch have an indentation in the central boundary which divides the rolls, while those 
on the unstable branch have a protrusion. Two-roll states have also been observed in the numerical simulations at 
Ra = 14 200 by Ma et al. 0] and at Ra = 37 500 by Leong [6]. The two-roll branches are very robust: both exist at 
least until Ra = 30 000 and the upper branch is stable for Ra < 28 086. We have been able to compute the CO branch 
only for 7167 < Ra < 10 348, and it is stable for Ra < 10087. The two-roll states have two symmetry axes (i.e. 
while the CO states, containing two light regions, one curved and one oval) have only one symmetry axis {Z2). 

Figure [19] shows that the high Rayleigh number {Ra > 20 000) asymmetric three-roll states and four-roll states 
greatly resemble two-roll states. Along all of the branches, the convective structures widen as Ra increases for all the 
branches: this is the form taken in this confined geometry of the well-known increase in wavelength for large systems 
of parallel rolls. For the branches emerging from the m = 1 and m = 2 primary bifurcations, this tendency eventually 
leads, after secondary bifurcations and more gradual deformations, to states which primarily contain two rolls. These 
are far removed from the trigonometric forms of the dipole and pizza states that prevail along these branches at low 
Ra. 



V. CONCLUSION 

We have presented an intricate bifurcation diagram describing Rayleigh-Benard convection in a cylinder with aspect 
ratio r = 2 and Pr = 6.7 for Ra < 30 000. This study is complementary to the time-dependent simulations described 
in our companion paper [4]; the branches of the bifurcation diagram were obtained by continuation from the stable 
states resulting from time integration. We have determined the bifurcation-theoretic origin of these states, including 
the five states observed experimentally by Hof et al. In one case, the path is straightforward: the four-roll branch 
results from a secondary bifurcation from the pizza branch, which in turn arises from a primary m = 2 bifurcation 
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FIG. 17: (Color online) Partial bifurcation diagram containing the two-roll and CO branches. Two-roll branches (solid, blue) 
arise via saddle-node bifurcation at Ra = 8677. CO branch (dashed, red) has been computed for 7167 < Ra < 10 348. The 
apparent intersections (between the two-roll and CO branches and between the stable and unstable two-roll branches) is an 
artifact of the projection. 
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FIG. 18: (Color online) Schematic partial bifurcation diagram showing the two-roll and CO branches. The two-roll branches 
originate at a turning point at Ra = 8677 and have been computed for Ra < 30 000. The upper branch is stable for Ra < 28 086. 
The CO branch has been calculated in the range 7167 < Ra < 10 348; it is stable for Ra < 10 087. 



from the conductive branch. In another case, it is more tortuous: for m = 3, two additional saddle-node bifurcations 
must be traversed between the primary marigold branch and the stable Mercedes branch which is actually observed. 
The torus and two-roll branches turned out to be disconnected (as far as we can tell) to the conductive state. Finally, 
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FIG. 19: (Color online) Four-roll (dashed, turquoise), asymmetric three-roll (solid, brick) and two-roll branches (solid, blue) at 
high Ra. 



we have located the primary m = 1 bifurcation leading to the three-roll branch, but it is accompanied by another 
simultaneously bifurcating branch and has unexpected stability properties. We have also traced the disconnected CO 
branch and the two-torus branch arising from a primary m = bifurcation. A schematic version of the bifurcation 
diagram is given in figure [20l while tables |TT] and |TTI1 list all of the branches we have obtained, as well as the bifurcations 
and their nature. 

The diagram we have obtained contains 17 branches of steady states, but is nonetheless incomplete. Although 
we have followed the primary branches originating at 5 bifurcations along the conductive branch, there are literally 
hundreds of other primary bifurcation points in the range Ra < 30 000. Each of the primary branches thus engendered 
can and does undergo many secondary bifurcations. In addition, while calculating the stability of the various branches 
we have observed many eigenvalues cross zero, signalling the appearance of a new branch. Finally, there is no way 
to ascertain how many other disconnected branches. It is surely unfeasible and unproductive to strive to find all 
branches. 

Despite the complexity of the bifurcation diagram, its main features can be described quite simply. Circle pitchfork 
bifurcations to trigonometric branches dictated by the geometry - dipole, pizza, marigold, two-torus - take place at 
low Rayleigh numbers. These undergo various secondary bifurcations at intermediate Rayleigh numbers that lead to 
states with rolls. At high Rayleigh numbers, there are three types of stable branches: torus, Mercedes, and states 
essentially containing two rolls. 

Our goal of relating the states obtained by time integration to the bifurcation diagram is largely realized, but there 
remain a few loose ends. One is to complete our understanding of the two simultaneously bifurcating m = 1 branches. 
Another small and clear-cut goal is to incorporate the time-periodic rotating S state we have described in jsl, Q into 
the bifurcation diagram by ascertaining its bifurcation-theoretic origin and exact domain of existence and stability. 
Comparing our study to that of Ma et al., we find very similar values for the five primary pitchfork bifurcations, for 
the secondary pitchfork bifurcation creating the four-roll branch, and for the two secondary pitchfork bifurcations 
which stabilize the two-torus branch. The other bifurcations in Table ITTTl are not present in Ma et al. These authors 
did, however, compute an additional stable four-spoked pattern at Ra = 14 200, whose bifurcation-theoretic genesis 
could be interesting to determine. 

Our study also points in several larger directions. It would be desirable to incorporate improved versions of the 
numerical methods we have used, namely adjustment of Ra within the Newton iteration, the inverse Arnoldi method, 
and the calculation of traveling waves as steady states in a rotating frame. This example could also be used as 
a test case to try to understand and to control the enormous variability in performance of BiCGSTAB in solving 
the linear equations of Newton's method in different regions of the bifurcation diagram. Finally, an extensive but 
straightforward goal is to compute a bifurcation diagram for the case in which the sidewalls are thermally conducting 
rather than insulating. Time-dependent simulations have already yielded initial conditions and approximate stability 
ranges for many branches [3|, For the conducting case, as for the insulating case, construction of a complete 
bifurcation diagram is doubtless impossible, but the stable steady and time-periodic states could all be traced back 
to their bifurcation-theoretic origin from the conductive state as we have done here. 
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FIG. 20: Schematic bifurcation diagram. Arbitrary quantity of the vertical axis chosen to ehminate all but one intersection. 
Bold lines indicate stable portions of branches. 



The complexity of the bifurcation diagram we have computed is interesting in light of the recent computational 
discovery of large numbers of unstable solutions of wall-bounded shear flows, e.g. [16|. It is hypothesized that weak 
turbulence can be understood as chaotic trajectories, e.g. [l^, that visit in turn the vicinities of the various unstable 
branches, which are created at saddle- node bifurcations. Our study contributes two observations to this line of 
research. First, this example provides a reminder that the existence of a large number of unstable solutions is a 
typical property of the hydrodynamic equations. Second, our study underlines the fact that such multiplicity can 
occur in the absence of complicated dynamics. 

This study showcases our numerical methods for carrying out time-integration, branch continuation and linear sta- 
bility analysis by using a single code with several different high-level drivers. Newton's method has many advantages: 
it is much faster and much more precise, and can compute unstable states. Time integration remains, nevertheless. 
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absolutely essential for generating initial states, especially since several important branches are disconnected from 
the conductive state. Although our cylindrical Rayleigh-Benard computation is quite specific, it demonstrates what 
can be accomplished for three-dimensional nonlinear problems by combining matrix-free preconditioned numerical 
methods with dynamical systems theory. 









nin 


-R^max 


family 


branch 


existence 


Stability 


Stability 


existence 


2 


pizza 


1849.4 


1879 


2353 


19450 




four-roll 


2353 


2353 


22 660 


23130 





upper two-torus 


1861.6 


2300 


5438 


12711 




(Cdj) 

lower two-torus 


2328.0 






12711 




upper one-torus 


3076.4 


4918 


> 30 000 


> 30 000 




lower one-torus ^ 


3076.4 


— 




> 17857 


3 


1 ] 

marigold 


1985.3 


— 




> 19 695 




mitsubishi 


4103 


— 




18 762 




cloverleaf 


4634.2 


— 




io /oz 




mercedes 


4634.2 


5503 


> 30 000 


> 30 000 


1 


tiger 


1828.4 


— 





> 7936 




(11) 

three-roll 


1828.4 


3762 


20 393 


> 30 000 




upper asymmetric three-roll 


21077.7 


21078 


> 30 000 


> 30 000 




do 

lower asymmetric three-roll 


21 077.7 






22125 


other 


upper two-roll 


8677.0 


8677 


28 086 


> 30 000 




(1) 

lower two-roll 


8677.0 






> 30 000 




CO 


< 7167 


< 7167 


10 087 


> 10 348 



TABLE II: List of all convective branches computed, with lower and upper limits of existence and stability. Stable patterns 
(between limits in boldface) should be observable in experiments or time-dependent simulations. Inequalities (> or <) indicate 
a lower or upper bound for the corresponding Ra: calculation of the branch either terminated for an unknown reason, or 
continued past 30 000. Dashes ( ) indicate the lack of stability limits for branches which are unstable throughout. 
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family 


bifurcation 


Ra 


comments 


2 


circle pitchfork 


1849.4 


creates pizza branch 




pitchfork 


2353 


creates four-roll branch 




eigenvalue crossing 


1879 


stabilizes pizza branch 




turning point 


19 450 


terminates pizza branch 




eigenvalue crossing 


22 660 


destabilizes four-roll branch 




turning point 


23 060 


terminates four-roll branch 





pitchfork 


1861.6 


creates upper two-tori branch 




pitchfork 


2328.0 


creates lower two-tori branch 




turning point 


12 711 


terminates upper and lower two-tori branches 




eigenvalue crossing 


2116 


stabilizes upper two-tori branch against m = 2 eigenvector 




eigenvalue crossing 


2300 


stabilizes upper two-tori branch against m = 1 eigenvector 




eigenvalue crossing 


5438 


destabilizes upper two-tori branch against m = 1 eigenvector 





turning point 


3076.4 


creates upper and lower one-torus branches 




eigenvalue crossing 


3330 


stabilizes upper one-torus branch against m = 6 eigenvector 




eigenvalue crossing 


3438 


stabilizes upper one-torus branch against m = 2 eigenvector 




eigenvalue crossing 


4408 


stabilizes upper one-torus branch against m = 5 eigenvector 




eigenvalue crossing 


4582 


stabilizes upper one-torus branch against m = 3 eigenvector 




eigenvalue crossing 


4918 


stabilizes upper one-torus branch against m = 4 eigenvector 


3 


circle pitchfork 


1985.3 


creates marigold branch 




pitchfork 


4103 


creates mitsubishi branch 




turning point 


4634.2 


creates mercedes and cloverleaf branches 




turning point 


18 762 


terminates cloverleaf and mitsubishi branches 




eigenvalue crossing 


5503 


stabilizes mercedes branch 


1 


circle pitchfork 


1828.4 


creates tiger and three-roll branches 




eigenvalue crossing 


3762 


stabilizes three-roll branch 




eigenvalue crossing 


20 393 


destabilizes three-roll branch 




turning point 


21077.7 


creates asymmetric three-roll branch 




sub critical pitchfork 


22125 


terminates asymmetric three-roll branch 


other 


turning point 


8677.0 


creates two-roll branches 




eigenvalue crossing 


28 086 


destabilizes two-roll branch 




eigenvalue crossing 


10 087 


stabilizes CO branch 



TABLE III: List of all bifurcations located. Circle pitchfork bifurcation breaks axisymmetry, creating "circle" of new states 
parametrized by azimuthal phase. Pitchfork bifurcation breaks a reflection symmetry, creating two branches. Eigenvalue 
crossings are necessarily accompanied by bifurcations, whose nature we have not investigated. Comments such as "stabilizes", 
"destabilizes" , "creates" , "terminates" are to be interpreted in the direction of increasing Rayleigh number. 
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